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A Fokker-Planck equation for the Wigner function evolution in a noisy Kerr medium (x^^^ non- 
linearity) is presented. We numerically solved this equation taking a coherent state as an initial 
condition. The dissipation effects are discussed. We provide examples of quantum interference, 
sub-Planck phase space structures, and Gaussian versus non-Gaussian dynamical evolution of the 
state. The results also apply to the description of a nanomechanical resonator with an intrinsic 
Duffing nonlinearity. 



I. INTRODUCTION 

Nonlinear interaction of light in a medium provides a 
very useful framework to study various nonclassical prop- 
erties of quantum states of radiation. These nonclassical 
properties are usually associated with quantum interfer- 
ence and entanglement. The phase space Wigner distri- 
bution function description of quantum states of light is 
a powerful tool to investigate such nonclassical effects. 
With the help of the Wigner function one can simply 
visualize quantum interference. For example, a signa- 
ture of quantum interference is exhibited in the Wigner 
function by the non-positive values and sub-Planck struc- 
tures [l]. The non-positive Wigner function is a wit- 
ness of a nonclassicality and monitors a decoherence pro- 
cess of a quantum state, e.g. a photon-added coherent 
state in the photon-loss channel [2 , photon-subtracted 
squeezed state [3], Gaussian quadrature-entangled sin- 
gle photon subtracted light pulse [4 , two-photon Fock 
state i5j, odd photon number states superposition [6] 
and coherent states superpositions [7 . The Wigner func- 
tion has also been used for computing numerically the 
quantum-mechanical corrections to the classical dynam- 
ics of a nanomechanical resonator (its characteristic pat- 
tern, the interference fringes and negative values, served 
as a signature for a classical to quantum domain transi- 
tion) [8j ; it has been applied in the model of the dynamics 
of a nanoscale semiconductor laser [9 , to list only a few 
examples of its applications. 

The Kerr medium provides one of the simplest non- 
linearity for which there exists a simple analytic Wigner 
function p!Ol expression. The highly x^^"^ nonlinear sys- 
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tems generated a lot of interest recently due to their ap- 
plications e.g. to nondemolition measurements [11], quan- 
tum computing architectures [12 , single particle detec- 
tors [13 . The well known example of such a medium is an 
optical fiber. However, the nonlinearity is small in fibers 
and is often accompanied by other unwanted effects. En- 
hanced Kerr nonlinearity was studied in terms of electro- 
magnetically induced transparency [M] and was observed 
in Bose Einstein condensates [TF and cold atoms [16]. 
Recent proposals predict obtaining enormous Kerr non- 
linearity using the Purcell effect [17 , Rydberg atoms [13], 
and interaction of a cavity mode with atoms [18]. The 
first and last method is predicted to obtain the nonlinear- 
ity of nine orders of magnitude higher than natural self- 
Kerr interactions with negligible losses. The significant 
nonlinearity has also been observed for nanomechanical 
resonators [20 . 

Most of the investigations of the Wigner function of 
light have been made for steady state situations. The 
simplicity of the Kerr medium will allow to study the 
full time-dependent Wigner function dynamics with or 
without a quantum noise. 

In this paper we present a Fokker-Planck equation 
which determines the time evolution of the Wigner func- 
tion in a noisy x^^^ medium. We solve this equation 
numerically assuming a coherent state as an initial con- 
dition. We discuss first an ideal and then a dissipative 
Kerr medium. The results obtained for the ideal case 
reveal the quantum nature of the state under evolution. 
The coherent state, known as the most classical among all 
the pure states, becomes a non-Gaussian squeezed state 
after some time of interaction with the medium. For 
some specific times it becomes a finite superposition of 
other coherent states. An interference pattern with the 
negative values is clearly visible on the plots of its Wigner 
function. Due to the small value of the x^^^ nonlinearity 
and losses in optical fibers such phenomena have never 
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been observed for light; neither for any other system. 
However, it turns out that not ah quantum effects are 
washed out due to the decoherence. The Fokker-Planck 
equation ahows for a similar state evolution analysis for 
any other input state. 

The paper is organized as follows: in section |n| the 
Fokker-Planck equation for the Wigner function evolu- 
tion in a x^^^ medium is derived from the Master equa- 
tion obtained for a single mode of light density opera- 
tor. The equation is displayed in a polar coordinates. 
The decoherence effects: losses and thermal noise are 
included. The initial and boundary conditions for the 
Wigner function are set. In section [TTI| the numerical re- 
sults of the Wigner function evolution in a nondissipative 
^(3) medium are introduced and discussed. These results 
are obtained using three independent methods: comput- 
ing the Fokker-Planck equation and the other two equa- 
tions determining the Wigner function directly, which are 
obtained from its definition. Correspondence between 
the Wigner function negative values and zeros of the 
Q-function is noted. In section |IV| the influence of the 



decoherence process on quantum effects such as the in- 
terference pattern and the negative values of the Wigner 
function is analyzed. The technical limitations on the use 
of the Fokker-Planck equation are discussed. The Wigner 
function sub-Planck structure is shown in section |Vl In 
section IVH the numerical methods used in sections Hill and 
IV are compared and discussed. Finally some conclusions 
are presented. 



II. THE FOKKER-PLANCK EQUATION FOR A 
SELF-KERR INTERACTION 



The interaction Hamiltonian for a totally degenerate 
four-wave mixing process, e.g. in an optical fiber, is of 
second-order in creation and annihilation a light op- 
erators 



(1) 



where is a nonlinear constant proportional to x^^'^ [2T] . 
Please note that similar Hamiltonian describes a single 
nanomechanical resonator with and a being raising 
and lowering operators related to its position and mo- 
mentum operators, and proportional to the Duffing 
nonlinearity [19 . 

In a general case, including damping and thermal noise 
in the medium, a one-mode density operator evolution, 
both for light and for a nanomechanical resonator, is de- 
termined by the Master equation [22] 



dtPir) 



■^[H,p{t)] + ^ (Mr), at] + [a,p(T)at]) 



TN[[a,p{T)],a% 



(2) 



where r = — A^t is a unitless evolution parameter, t is 
the interaction time, F is a damping constant, N = 



l/{exp (1^) — l}isa mean number of photons in a ther- 
mal reservoir. 

The solutions of the above equation are well known [2TJ 
l23l [24[ |25] but since it is an operator- valued equation, 
they are inappropriate for computer simulations. A 
Fokker-Planck type evolution equation for the Wigner 
function can be easily obtained from equation ([2| using 
the standard quantum optics phase space methods [26 . 
Since every density operator determines its Wigner func- 
tion uniquely, the knowledge of its evolution is equivalent 
to the full knowledge of the density operator dynamics. 

The dynamics of the Wigner function in a dissipative 
medium with a self-Kerr interaction is governed by the 
Fokker-Planck equation, which takes the following form 
in the polar coordinates 7 = re*"^ 



drW{r,r,^) = Ur^-l)d^ 



1 /I 



1 



WiT,r,ip), (3) 



where 7 is a point in a phase space, ^ = T/k,. 

This is a third-order nonlinear differential equation. 
We compute this equation for the following initial and 
boundary conditions. The evolution starts with a coher- 
ent state \a) described by a Wigner function W(0,7) = 
lg-2|Q!-7| rpj^^ Wigner function tends to zero in the 
infinity lim^^^o W{r^ r, (f) = 0. Since some of the coef- 
ficients in equation (|3| are singular for 7 = 0, we take 
iy(r, 0,0) = |:e~^l"l ^ ^, which was derived from the 
Master equation. 



III. THE NONLINEAR PART OF EVOLUTION 

For ^ = and = the Fokker-Planck equation (|3| 
reduces to its first line and describes the evolution in an 
ideal Kerr medium. 

Figures [l]|2j|3j [4] present the plots of the Wigner func- 
tion for different parameters r and a = 5. Due to the fact 
that the number operator a^a is a constant of motion, the 
probability distribution during its evolution will remain 
situated within a circle of radius \a\. Beginning from a 
circular shape, the Wigner function turns into an ellipse 
and squeezing appears in an appropriate direction. Then 
the ellipse changes into a banana shape and a "tail" of the 
interference fringes appears where the distribution takes 
the negative values. Squeezing increases and the state 
becomes non-Gaussian. The "banana part" of the distri- 
bution gets thinner in the radial direction and becomes 
more and more smeared in the azimuthal direction. 

The rotation and stretching effects which lead to 
squeezing are revealed by the first term in the first line of 
the equation (|3| . It shows the dependence of the angular 
velocity on the distance r from the origin of coordinates. 
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This corresponds to the fact that the nonhnear refrac- 
tive index 72(7) of Kerr medium is intensity dependent 
77,(7) = ^0 + ^2|7p- The intensity fluctuations modulate 
the nonlinear refractive index and this in turn modulates 
the phase of traveling light. Photons with stronger am- 
plitude will acquire phase faster than the photons with 
smaller amplitude. 

The other first line terms of this equation consisting of 
mixed and third order derivatives are responsible for the 
interference fringes formulation. 

If r is taken as a fraction of the period of the evolution, 
r = 27tR where < 1 is a rational number, the initial 
coherent state becomes a superposition of other coherent 
states, of the same amplitude but different phases, see 
fig. [4] [27 . This is is also known as a fractional revival [28j . 

The evolution governed by the self-Kerr interaction 
Hamiltonian is periodic. It is easy to see this using a 
unitary evolution operator U{r) = exj[){—iHt/h} rather 
than the Fokker-Planck equation approach. Its action on 
a coherent state \a) expressed in the Fock state basis 



|^(a,r)) = /7(r)|a) = e-^y ^ e^^^^'-^V) (4) 



n=0 



shows that |^(a,r)) = |^(a,r + 27r)). In particular, for 
r = 2/c7r where k is an integer, equation Q reduces again 
to the coherent state \a). The same periodicity holds for 
the Wigner function. This condition has been used as an 
additional check point in our numerical simulations. 

For the wave function given by equation Q one can 
determine the Wigner function as an analytical function 
of the evolution parameter r using the definition [29] 



W{t, 7) = ^e^^-^^" j<fp{-(3\p{a, t)W 



-2(/37*-r7) 



(5) 

where p{a,r) is the density operator of the state 
r)). We apply the above equation and express it in 
two equivalent ways which have been used to obtain the 
numerical results for the nondissipative Kerr medium 



W(r,7) = -e-^l^l e 

TT 



e "2 



^ (2a7*e- 

k=0 



q=0 



ql 



k\ 



(6) 



n,m=0 

(7) 

The equations Q and ^ explain how the interference 
fringes arise from the interference of the Fock number 
states \n). By virtue of equation ([7|), the Wigner function 
can be viewed as a sum of the two components: the first 
one, independent of r, being a superposition of the Fock 



states' Wigner functions (7) and the second, being 

an evolution dependent component /(r, 7), responsible 
for the evolution of the fringes 



E^W^HH(7) + /(r,7). (8) 



n=0 



Although the equation (|6| and Q are special solutions 
of the Fokker-Planck equation, we obtained them sepa- 
rately to proceed independent computation. 

The effects of squeezing and error contour rotation 
in the phase space can also be observed by study- 
ing the easy-to-compute Q-function evolution (5(r, 7) = 
l/7r(7|/)(a,r)|7) 



Q(r,7) 



1 



-I«|^-|7l^ 



n 



(9) 



The contour plot of the Q-function is plotted in ^ . So- 
lutions of the Fokker-Planck equation for the Q-function 
in a dissipative and noisy Kerr medium have been 
widely studied EH |23l [H [25]. Negativities achieved 
in the Wigner function correspond to zeros of the Q- 
function [30 . 

All the figures [l]-[4] have been obtained computing the 
equations ([3|, ([6| and ^ independently. 



IV. DISSIPATION EFFECTS 

In the presence of damping the Wigner function will 
not remain situated within at the circle of the radius \a\ 
any more. The second line in the Fokker-Planck equa- 
tion, proportional to the first radial derivative, makes 
the Wigner function moving towards the origin of the co- 
ordinates. This effect corresponds to decreasing number 
of quanta in the state during the evolution. At the origin 
the state becomes vacuum. 

The third line describes the effects of dispersion. 

It is interesting to note an interplay between the non- 
linear evolution and damping terms in the Fokker-Planck 
equation. If damping is negligible, the numerical simu- 
lations obtained for small \a\ :^ 1 and large |a| 1 
initial amplitudes are very similar. At the beginning of 
evolution the squeezing effect dominates and then the 
interference appears. It is not the case if the losses are 
included in the simulations. If the amplitude is too small 
the decoherence washes out all the quantum effects and 
the Gaussian Wigner function located at the coordinates' 
origin, genuine to vacuum state, is obtained very quickly. 

We proceed the numerical computation for experimen- 
tally reasonable values genuine to a nanomechanical res- 
onator. The nonlinear constant n c^2 • 10~^s~^ and the 
damping constant F :^ 10^ s~^ = 10^). The ther- 
mal noise coefficient in the room temperature is equal to 

The nonlinear evolution term for r :^ becomes 
meaningful \i r \a\. li \a\ 1 the damping term 
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^ :^ 0.5 • 10^ dominates over the nonhnear :^ 1 and 
interference terms. The coherent state wih be almost 
immediately radially displaced towards the origin of co- 
ordinates. For \a\ 0^ 10^ all the effects are in balance: 
r^c^lO^o and §^ c^O.S-lO^o. 

Computing the Fokker-Planck equation for such large 
values of coefficients {a 10^ ^ 10^ = 1.9 • 10-^"^) 
would require a very dense and large grid, thus a lot of 
computer memory. Therefore, to visualize an influence 
of these two effects on the Wigner function, we rescaled 
the parameters. 

As we pointed out above, in order to see the nonlinear 
effects during the Wigner function evolution in presence 
of dissipation, the input of nonlinear term has to be of 
the same order as the dissipation term is. Therefore, 
we keep the ratio between them constant and equal to 
the ratio evaluated for the non-rescaled case: (, ^ a and 
N = 1.9 ■ For a = 2 they are equal to ^ 

2 and N = 3.8 • 10~^^. Although such a small value 
of the thermal coefficient implies a negligible effect on 
the evolution, we included it in the simulation. Such 
simulations are accessible using a standard PC computer 
with 4 GB of memory. 

Figures [Sj [6) [7| [S] present the numerical simulations of 
the Fokker-Planck equation obtained for a = 2, A/" = 3.8- 
10~^^ and the following values of the damping constant: 
^ = 0, ^ = 0.1, C = 1,^ = 2. 

Due to the fact that damping washes out the inter- 
ference effects, any nonzero damping constant prevents 
both: coherent superposition state formulation and pe- 
riodicity of the evolution. As it is well known [31 these 
states are extremely fragile to the decoherence. The ar- 
eas where the Wimer function takes the negative values 
disappear, see fig. [sj For an ideal x*^^^ medium the neg- 
ativities show up for 0.08 + 2/c7r < r < 6.20 + 2/c7r where 
k is an integer. If ^ = 0.1 the negativities are present 
for 0.08 < r < 5.80. However, they do not appear in the 
second round (0.08 + 27r < r < 5.80 + 27r). If ^ = 1 the 
negative values are present only for 0.08 < r < 1.52 and 
for 0.08 < r < 0.89 if ^ = 2. 

The Wigner function obtained for the coherent super- 
position states evolution times r = , r = tt and r = 27r 
is depicted on figures [6j [7| |8] For ^ = 0.1 the structure of 
superposition coherent states is still well preserved. How- 
ever, some additional circular trails are present. These 
are specially visible for r = tt and r = 27r. For ^ = 2 and 
T = ^ the state is already a vacuum. The vacuum state 
is achieved for r = tt (r = IOtt) if ^ = 1 = 0.1). 



V. SUB-PLANCK STRUCTURE IN 
PHASE-SPACE 

In the framework of quantum mechanics, it has been 
recognized that small structures on the sub-Planck scale 
do show up in quantum linear superpositions [32l l33] . 
Such sub-Planck structures can be shown, if the Wigner 
function is plotted in a phase-space region below the 



Heisenberg relation. That means in the phase space 
volume less than 1, since AX1AX2 > 1 for the ampli- 
tude Xi and phase X2 quadratures. The quadratures 
are given by the annihilation a = ^{Xi + iX2) and cre- 
ation = ^{Xi — 1X2) operators and their mean values 
correspond to phase-space coordinates in the following 
way (Xi) = Re 7, {X2) = Im7. 

Figure [9] compares behavior of the sub-Planck struc- 
ture of the Wigner function obtained for the compass 
state (r = 7r/2) and the Schrodinger cat state (r = tt) 
for a = 2 in the ideal x^^^ medium and in presence of 
damping = 0.1). The regions of the opposite sign, 
in the form of "dots," are clearly visible. They become 
flattened and smeared by the damping effects. 

The sub-Planck structure is also present for the ratio- 
nal values of evolution parameter, for which the state Q 
does not correspond to any coherent superposition state. 
It is depicted on figure 10 for r = 0.16, r = 0.3, r = 0.6, 
r = 1 and a = 5. This structure is of different form. The 
areas of positive and negative values form "ribbons." 



VI. COMPARISON OF THE NUMERICAL 
METHODS 



Figures [T]-[4] presenting the Wigner function evolution 
in an ideal x^^^ medium, have been obtained by three 
independent numerical simulations: computing of the 
equations ([3|, (|6| and ([7). The numerical results of the 
evolution in a noisy medium have been achieved using 
only the Fokker-Planck equation 

All programs used to simulate the Wigner function 
were written in standard C++ language, due to its high 
speed and good portability. The structure of the pro- 
grams used for computation of the equations (|6| and ^ 
was based strictly on these equations: they both consist 
of two loops which iterate on q and k (or m and n, re- 
spectively) and sum up the evaluated coefficients for a 
given r and 7. Common terms were moved out from the 
loops (especially the inner one) to optimize the evalu- 
ation of the sums. The main function was enclosed in 
input /output logic which calls it for different values of 7 
and writes the results to a file. 

The equations ([6| and ^ involve summation of the 
infinite series. However, only a finite number of their 
terms can be computed. Therefore, assuming a given 
precision of the simulation, we cut the series off. 

The series in equation (|6| is much slower convergent 
than the series in ^ because it contains fast oscillating 

exponential term e~l"l^^'^^^ Depending on the phase 
factor in the exponent, this term can take either large or 
small values, which makes impossible to compute equa- 
tion ([6| using the standard double precision floating point 
number representation. Therefore, for \a\ < 5 both infi- 
nite sums, internal and external, are cut off by including 
at least 500 terms in each of them. Also very high pre- 
cision (at least 25 significant digits) needs to be applied 
for the computation. Such a high precision was obtained 
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using the free Class Library for Numbers (CLN) and its 
long float type. 

On contrary, it is enough to use standard double preci- 
sion numbers to compute equation However, we also 
applied the CLN here due to the fact that it can deal 
with the complex numbers while C++ itself cannot. The 
same precision was used while computing equation ^ 
however, it was enough to take 100 terms in both sums 
into account to obtain the same results as these which 
were achieved by computing This made simulations 
of equation ([7| much faster. Additionally, the product of 
the exponens derivatives was replaced with the following 
formula 34 



k=0 



-4i7r 



i(n,m) 



-4i7r 



E 

fc=0 



n\ ml 



kj {m — k)\ 



x(-47)^-^(-47*)^ 



(10) 



Both methods, equations (|6| and ([7|), allow for compu- 
tation of the Wigner function for an arbitrary evolution 
parameter r and point 7 in the phase space. Evaluation 
of the function for a chosen r does not require simu- 
lation of the whole evolution. To visualize the Wigner 
function we calculated it for each point of a 100 x 100 
square mesh of the size proportional to the value of a 
(—5 ... 5 for a = 2, —8 ... 8 for a = 5, — ^ • • • ^ for sub- 
Planck structures), to show all interesting aspects of the 
simulated function. 

On the contrary, numerical simulation of the Fokker- 
Planck equation (|3| involves solving a partial differential 
equation (PDE), which in turn requires setting a grid 
(discretized phase space) and solving a set of linear equa- 
tions (one equation for each point of the grid) at each 
time step. 

To solve the PDE, the equation (|3| was rewritten to 
difference quotient form, where all the partial deriva- 
tives were replaced with appropriate difference quo- 
tients of fourth order to keep the order of the quo- 
tients even and higher than the order of the equation 
at the same time. Even difference quotients are more 
convenient, since they are symmetric with respect to 
the point where the derivative is taken. The quotients 
were determined using the Mathematica program and its 
FiniteDif f erencePolynomial function. We used polar 
mesh coordinates with constant radial distance Ar and 
constant angle A(p between points. Below we present the 
first, second and third derivatives discretized up to the 
forth order: 



dx- 



fk-2 - 8/fc-i + 8/fe+i - /fc+2 
12h 



d: 



,2 ^ -fk-2 + 16/fc-l - 30/fe + 16/fe+i - /fc+2 

12/^2 

-fk-2 + 2/fc-i — 2/fc+i + /fc+2 

2h^ 



(11) 

(12) 
(13) 



where x = r^(f, h = Ar, A(/? and k = i,j for r and (f 
respectively, fk is the Wigner function value. The mixed 
derivatives where built up combining the basic above for- 
mulas. 

Since the standard explicite computation methods 
turned out to be unstable due to large number of time 
steps, we applied implicite method. At the same time, 
it also allowed for larger time steps. The main idea 
of implicite method is to determine the time derivative 
between steps r and r + Ar using the spatial deriva- 
tives obtained from the next time step W(r + Ar, 7, 7*), 
not from the previous one VK(r, 7,7*). This method 
requires solving a set of linear equations, with Wigner 
function values in the next time step W{r + Ar, 7,7*) 
for each point 7 of the mesh as unknowns. The 
stability of this method relies on well-chosen ratio of 
the time step and the grid density. For a = 2 we 
used time step Ar = 3^ and the grid consisting of 
300 (radial direction) x 540 (azimuthal direction) points. 
Thus, there were 162 • 10^ points and equations in the 
set. 

In order to solve such a large set of equations, we put 
it in a matrix form. To achieve it we assigned an index 
i for each point of the polar mesh using the formula i = 
Nfr/Ar + cp/Acp^ where Nf is a density of the grid in 
azimuthal direction, r and cp are polar coordinates of a 
point of the mesh and Ar and A(f are radial and angular 
distances between mesh points. This index was then used 
as a coordinate to the row and column of the matrix 
where the coefficients of the equations were placed and 
to rewrite the Wigner function values of the previous time 
step from grid to vector form and the result from vector 
back to mesh form. 

The index i (described above) was chosen in such a 
way that the matrix of equation coefficients is a sparse 
band matrix. The sample of such matrix computed for 
small mesh size is depicted on fig. [TT] Black points rep- 
resent nonzero items (single numbers). It consists of five 
diagonal five-element bands and additional two includ- 
ing either two or one element. In practice, the matrix is 
much larger and has 162 • 10^ x 162 • 10^ elements, but the 
characteristic band diagonal structure with five-element 
bands is preserved. 

The set of equations is solved using Band Diagonal 
Systems method [35]: the matrix is inverted applying 
LU decomposition (Crout method) and then solved with 
Gaussian elimination while keeping the matrices in spe- 
cial compressed form to save computer memory. We used 
the bandec and banmul routines from [3F, translated to 
C++ and optimized. We also used standard double pre- 
cision floating point numbers; there was no need for CLN 
library. The main computing routine was accompanied 
by input /output logic used to read simulation parameters 
and write out the results. 

The flow of the program consisted of an introductory 
step - performing a LU decomposition and calculating 
Wigner function for r = (initial values for each point of 
the grid) - and a sequence of steps calculating evolution 
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of the Wigner function for each r = tq + nAr. The 
solution gives the values of the function for all 7 at each 
evolution step. 

During the simulation some sanity checks were taken 
to make sure that the calculations are done properly, e.g. 
the integral of the Wigner function over the phase space 
is equal to 1 in each time step. The value of the integral 
unequal to 1 was a sign of either too small mesh density 
or too big time step or too small precision of the calcu- 
lations. We also compared the results of simulation of 
equations ([s]), ([6| and ^ for the same input parameters. 

The numerical results are presented using OpenDX 
program. 

VII. CONCLUSION 

In this paper we presented the Fokker-Planck equa- 
tion which allowed for numerical computation of the full 



Wigner function evolution governed by the self-Kerr in- 
teraction Hamiltonian. This equation can be used for a 
state analysis for any input state, mixed or pure. We 
took a coherent state as an input state for the simula- 
tion. For a decoherence-free evolution the results have 
been obtained by three different numerical algorithms. 
We discussed the influence of the decoherence process 
on the nonclassicality of the state under evolution. For 
an exemplary calculation we took into consideration the 
experimentally reasonable values of the nonlinearity and 
damping constant for a nanomechanical resonator. We 
also presented the sub-Planck structure of the Wigner 
function which arises during the evolution. This equation 
can be applied to any system described by the self-Kerr 
interaction. 
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FIG. 1: (Color online) The Wigner function for a = 5 and r = is a Gaussian function for the coherent state - the left figure. 
Due to the fact that the evolution in an ideal x^^^ medium is periodic, we achieve the same shape of the Wigner function for 
r = 27T. The distribution is an ellipse and the state becomes squeezed beginning from r = 0.01 - the middle figure. The 
negative values of the function start to appear for r = 0.04 - the right figure. 
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FIG. 2: (Color online) Further steps of Wigner function evolution obtained also for a = 5 and r = 0.08 - the left figure, 
T = 0.16 - the middle figure and r — 0.3 - the right figure. An interference pattern arises in the form of a "tail" of interference 
fringes. 




FIG. 3: (Color online) For the evolution parameter values close to r = 27ri?, where < 1 is a rational number, in some ares a 
destructive and in the other a constructive interference starts to dominate - the left (r = 0.6) and the middle figure (r = 1), 
thus leading to a Wigner function of a coherent superposition state - the right figure (r = 7v/3). Figures evaluated for a = 5. 
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FIG. 4: (Color online) The superposition states obtained for a = 5 and r = 27r/5 - the left figure, r = 7r/2 (the compass state) 
- the left middle figure, r = 27r/3 - the right middle figure and r = tt (the Schrodinger cat state) - the right figure. 
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FIG. 5: (Color online) The effect of the decoherence on the Wigner function evaluated for a — 2, r = 0.27r, the thermal noise 
N — 3.8 • 10~^^ and the different values of the damping constant: ^ = - the left figure, <J = 0.1 - the left middle figure, ^ = 1 
- the right middle figure, ^ = 2 - the right figure. Note that the thermal noise gives no input if ^ = 0. 
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FIG. 6: (Color online) The coherent superposition states are extremely fragile to the decoherence process. The Wigner function 
of four coherent state superposition (r = 7r/2) of the amplitude a — 2 evaluated including the thermal noise N — 3.8 • 10~^^ 
and the damping constant equal to ^ = - the left figure, ^ = 0.1 - the left middle figure, <f = 1 - the right middle figure, 
^ = 2 - the right figure. 
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FIG. 7: (Color online) The decoherence looks similarly for the two coherent state superposition Wigner function (r = tt), also 
of amplitude a = 2, the same thermal noise N = 3.8 • 10~^^ and the damping constants: ^ = - the left figure, ^ = 0.1 - the 
left middle figure, <J = 1 - the right middle figure, ^ = 2 - the right figure. 
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FIG. 8: (Color online) The Wigner function evolution is a noisy x medium is not periodic any more. The distribution 
evaluated for a = 2, r = 27r, the thermal noise N = 3.8 • 10~^^ and the damping constant equal to ^ = - the left figure, 
^ = 0.1 - the middle figure, <J = 1 - the right figure. The state pursues vacuum. 




FIG. 9: (Color online) A sub-Planck structure of the Wigner function for r = 7v/2 and ^ = - the left figure, r = 7r/2 and 
^ = 0.1 - the left middle figure, r = tt and ^ = - the right middle figure, r = tt and ^ = 0.1 - the right figure. 
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FIG. 11: The sparse band matrix necessary to invert to solve the set of Unear equations giving the Wigner function values for 
all points of the discretized phase space. 



